#----this is the R version of the ML for the erator ----
require(irr) # for kappa2
require(rpart) # for CART
require(randomForest)
require(e1071) # for svm
require(kernlab) # for gausspr
require(neuralnet) # for neural network
require(kknn) # for k nearest neighbor

#---defining some evaluation functions -----

roundScore<-function(escore){
  escore[escore<1.5] = 1
  escore[escore>=1.5 & escore<2.5]=2
  escore[escore >= 2.5 & escore < 3.5] = 3
  escore[escore >= 3.5 & escore < 4.5] = 4
  escore[escore>= 4.5 & escore < 5.5] = 5
  escore[escore>= 5.5] = 6
  return(escore)
}

PercentageExactMatch <-function(trueScore,fittedScore){
  nmatch = length(trueScore[trueScore == fittedScore])
  ntotal = length(trueScore)
  return(nmatch/ntotal*100) 
}

standDiff<-function(h1,escore){
  mh1 = mean(hscore)
  me = mean(escore)
  res = (mh1 - me)/sqrt((var(hscore)+var(escore))/2)
  return(res)
}

evaluateResult<-function(hscore,escore){
  escoreRND = roundScore(escore)
  print(" --- kappa ---")
  print(kappa2(cbind(hscore,escoreRND),weight="squared")$value)
  print("---Percentage of exact match ---")
  print(PercentageExactMatch(hscore,escoreRND))
  print("-----cross match table ----")
  print(table(hscore,escoreRND))
  print("--- mean escore ---")
  print(mean(escoreRND))
  print(" --- std escore ---")
  print(sd(escoreRND))
  print("--- standardized difference ---")
  print(standDiff(hscore,escoreRND))
  print("----correlation hscore and escore--")
  print(cor(hscore,escore))
}


# --- read in training data ---
#tdata = read.table('E:/rapidminerData/MachineLearning/EraterResults/arggen10_mb.sas7bdat.csv',sep=',',header=T)
#tdata = read.table('/media/jghao/data/rapidminerData/MachineLearning/EraterResults/arggen10_mb.sas7bdat.csv',sep=',',header=T)

#tdata = tdata[tdata$h1 != 0,]
#tidx = sample(dim(tdata)[1],10000)
#tdata = tdata[tidx,]

#vdata = read.table('E:/rapidminerData/MachineLearning/EraterResults/arggen10_xv.sas7bdat.csv',sep=',',header=T)
#vdata = read.table('/media/jghao/data/rapidminerData/MachineLearning/EraterResults/arggen10_xv.sas7bdat.csv',sep=',',header=T)

tdata = read.table('/home/jghao/research/LPworknew/06222013/arg_mb.csv',header=T,sep=',')
vdata = read.table('/home/jghao/research/LPworknew/06222013/arg_xv_noflag.csv',header=T,sep=',')
tdata = tdata[tdata$h1 != 0,]
vdata = vdata[vdata$h1 != 0,]
hscore = vdata$h1

#----------------------
#---Linear model ------
#----------------------

m_lin = lm(h1~LOGDTU+LOGDTA+nsqg+nsqm+nsqu+nsqStyle+nWF_MEDIAN+WORDLN_2+ColprepSVF+DWU,data=tdata)
escore = predict(m_lin,vdata)

#---evaluate linear model ----
evaluateResult(hscore,escore)


#----------------------
#---CART regression ---
#----------------------

m_cart = rpart(as.character(h1)~LOGDTU+LOGDTA+nsqg+nsqm+nsqu+nsqStyle+nWF_MEDIAN+WORDLN_2+ColprepSVF+DWU,data=tdata)
escore = predict(m_cart,vdata)

#---evaluate CART (regression) ----
evaluateResult(hscore,escore)
cor(vdata$ARGSCORE,escore)
plot(vdata$ARGSCORE,escore,xlab='HScore',ylab='EScore',pch=',',main='CART Regression')

#----------------------
#-CART calssification -
#----------------------

m_cart = rpart(h1~LOGDTU+LOGDTA+nsqg+nsqm+nsqu+nsqStyle+nWF_MEDIAN+WORDLN_2+ColprepSVF+DWU, method='class',data=tdata)
escore = predict(m_cart,vdata)

#---evaluate CART (calssification) ----
evaluateResult(hscore,escore)

#------------------------
#--k nearest neighbor----
#------------------------

fit = kknn(h1~LOGDTU+LOGDTA+nsqg+nsqm+nsqu+nsqStyle+nWF_MEDIAN+WORDLN_2+ColprepSVF+DWU,tdata,vdata, k=20, kernel="rectangular")
escore=fitted(fit)
evaluateResult(hscore,escore)
cor(vdata$h1,escore)


#plot(vdata$ARGSCORE,escore,xlab='HScore',ylab='EScore',pch=',',main='Random Forest')

require(FNN)
res=knn.reg(train=tdata[,3:12],test=vdata[,3:12],tdata$h1,k=20,algorithm="kd_tree")
escore=res$pred
evaluateResult(hscore,escore)
cor(vdata$h1,escore)

#------------------------
#--random forest --------
#------------------------

fit = randomForest(h1~LOGDTU+LOGDTA+nsqg+nsqm+nsqu+nsqStyle+nWF_MEDIAN+WORDLN_2+ColprepSVF+DWU,data=tdata,ntree=1200,mtry=6)
escore=predict(fit,vdata)
evaluateResult(hscore,escore)
cor(vdata$h1,escore)


#plot(vdata$ARGSCORE,escore,xlab='HScore',ylab='EScore',pch=',',main='Random Forest')


#------------------------
#--        SVM   --------
#------------------------

x = cbind(tdata$LOGDTU,tdata$LOGDTA,tdata$nsqg,tdata$nsqm,tdata$nsqu,tdata$nsqStyle,tdata$nWF_MEDIAN,tdata$WORDLN_2,tdata$ColprepSVF,tdata$DWU)
y = tdata$h1
vx = cbind(vdata$LOGDTU,vdata$LOGDTA,vdata$nsqg,vdata$nsqm,vdata$nsqu,vdata$nsqStyle,vdata$nWF_MEDIAN,vdata$WORDLN_2,vdata$ColprepSVF,vdata$DWU)


fit =  svm(x, y, type='eps-regression')
escore=predict(fit,vx)
print ('svm default')
evaluateResult(hscore,escore)

fit =  svm(x, y, degree = 2)
escore=predict(fit,vx)
print ('svm: 2 degree kernel ')
evaluateResult(hscore,escore)

fit =  svm(x, y, degree = 4)
escore=predict(fit,vx)
print ('svm: 4 degree kernel ')
evaluateResult(hscore,escore)

fit =  svm(x, y, degree = 8)
escore=predict(fit,vx)
print ('svm: 8 degree kernel ')
evaluateResult(hscore,escore)

fit =  svm(x, y, kernel = "polynomial",degree = 3)
escore=predict(fit,vx)
print ('svm: 3 degree polynomial kernel ')
evaluateResult(hscore,escore)

fit =  svm(x, y, kernel = "polynomial",degree = 5)
escore=predict(fit,vx)
print ('svm: 5 degree polynomial kernel ')
evaluateResult(hscore,escore)

#---best ---
fit =  svm(x, y, gamma = 0.02,cost = 15)
escore=predict(fit,vx)
print ('svm default')
evaluateResult(hscore,escore)
cor(vdata$h1,escore)


#cor(vdata$ARGSCORE,escore)
#plot(vdata$ARGSCORE,escore,xlab='HScore',ylab='EScore',pch=',',main='SVM')


#------------------------
#-- Gaussian process-----
#------------------------

x = cbind(tdata$LOGDTU,tdata$LOGDTA,tdata$nsqg,tdata$nsqm,tdata$nsqu,tdata$nsqStyle,tdata$nWF_MEDIAN,tdata$WORDLN_2,tdata$ColprepSVF,tdata$DWU)
y = tdata$h1
vx = cbind(vdata$LOGDTU,vdata$LOGDTA,vdata$nsqg,vdata$nsqm,vdata$nsqu,vdata$nsqStyle,vdata$nWF_MEDIAN,vdata$WORDLN_2,vdata$ColprepSVF,vdata$DWU)

fit =  gausspr(x, y)
escore=predict(fit,vx)

evaluateResult(hscore,escore)


#cor(vdata$ARGSCORE,escore)
#plot(vdata$ARGSCORE,escore,xlab='HScore',ylab='EScore',pch=',',main='Gaussian Process')


#------------------------
#-- Neural Network  -----
#------------------------

x = cbind(tdata$LOGDTU,tdata$LOGDTA,tdata$nsqg,tdata$nsqm,tdata$nsqu,tdata$nsqStyle,tdata$nWF_MEDIAN,tdata$WORDLN_2,tdata$ColprepSVF,tdata$DWU)
y = tdata$h1
vx = cbind(vdata$LOGDTU,vdata$LOGDTA,vdata$nsqg,vdata$nsqm,vdata$nsqu,vdata$nsqStyle,vdata$nWF_MEDIAN,vdata$WORDLN_2,vdata$ColprepSVF,vdata$DWU)

fit=neuralnet(h1~LOGDTU+LOGDTA+nsqg+nsqm+nsqu+nsqStyle+nWF_MEDIAN+WORDLN_2+ColprepSVF+DWU,data=tdata,hidden=10)

escore=predict(fit,vdata)

evaluateResult(hscore,escore)
cor(vdata$h1,escore)

plot(vdata$ARGSCORE,escore,xlab='HScore',ylab='EScore',pch=',',main='Gaussian Process')
